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ABSTRACT 


A new  method  for  the  numerical  simulation  of  three-dimensional 
incompressible  flows  is  described.  Our  vortex-in-cell  (VIC)  method 
traces  the  motion  of  the  vortex  filaments  in  the  velocity  field  these 
filaments  create  on  an  Eulerian  mesh  via  the  fast  integration  of  a 
Poisson  equation.  By  incorporating  the  viscous  or  subgrid-scale 
effects  into  a filtering  procedure,  the  computed  scales  of  motion 
are  assumed  to  be  essentially  inviscid.  Results  on  tracing  a 
periodic  array  of  single  vortex  rings  are  compared  with  a Green's 


function  calculation. 


1.  INTRODUCTION 


In  this  paper,  we  describe  a new  method  for  the  numerical  simula- 
tion of  three-dimensional  incompressible  flows.  Our  approach  differs 
from  other  numerical  fluid  simulations  in  that,  rather  than  solving 
the  Navier-Stokes  equations  on  an  Eulerian  mesh,  we  trace  the  motion 
of  the  vortex  filaments  in  the  velocity  field  these  filaments  create. 

In  addition,  the  velocity  field  is  not  calculated  directly  by  Biot- 
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Savart's  law  of  interaction  as  inL  * * , but  by  creating  a mesh-record 

of  the  vorticity  field,  then  integrating  a Poisson  equation  to  get  the 
stream  function  and  generating  a mesh-record  of  the  velocity  field. 

The  vortex  filaments  are  then  stepped  forward  in  time. 

Vortex  pushing  methods--as  distinct  from  Navier-Stokes  methods-- 
have  a history  of  success  in  two  dimensions^4  In  three  dimen- 

sions^’^’^  , one  of  us  has  applied  it  to  a small  number  of  simple 
vortex  filaments,  but  at  considerable  cost  because  of  the  time  required 
to  sum  all  the  mutual  Biot-Savart  type  interactions  between  the  many 
elements  in  all  the  filaments.  The  "cell"  or  "mesh"  method  speeds  up 
the  calculation  of  the  interactions  and  allows  the  three  d mensional 
vortex  pushing  method  to  be  applied  to  a space  densely  filled  with 
vortex  filaments,  each  filament  being  resolved  in  fine  detail  along  its 
length.  The  techniques  of  making  optimal  use  of  the  mesh  for  evaluating 
interactions  between  different  fluid  elements  are  the  same  as  those 

r9“12] 

used  by  plasma  physicists  for  calculating  particle  interactions^ 

In  creating  a code  for  the  evaluation  of  local  flow  fields  due  to  a 
family  of  vortex  filaments  we  have  therefore  taken  over  the  principles 
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and  architecture  of  a code  developed  for  magnetic  field  evaluations 

rial 

in  plasma  simulations1' 

An  outline  of  the  remaining  sections  of  the  paper  is  as  follows. 
The  fundamentals  of  vorticity  dynamics  relevant  to  our  technique  are 
discussed  in  Section  2 and  a description  of  the  computational  method 
is  given  in  Section  3«  1°  Section  4,  the  results  of  computational 

experiments  involving  vortex  rings  are  presented  and  discussed  while 
conclusions  and  suggestions  for  further  work  are  given  in  Section  5* 


2.  BASIC  PRINCIPLES 

We  assume  an  unbounded  incompressible  flow,  fully  periodic  in  the 
three  dimensions.  In  each  periodic  box,  the  vorticity  field  consists 
of  a collection  of  vortex  filaments.  The  governing  dynamic  equation 
for  the  vorticity  in  these  filaments,  w = V x u , is 


^ = Hi  • Vu  + VV^  U 


(2.1) 


where  v is  the  viscosity  and  the  velocity  field  is  determined  kinema- 
tically from 


u = - V X uj 


(2.2) 


In  the  following  we  assume  that  the  computed  scales  of  motion  are  essen- 
tially inviscid.  Any  viscous  or  subgrid-  scale  effects  on  the  computed 
fields  are  incorporated  into  the  filtering  procedure  introduced  below 
and  described  in  more  detail  in  the  next  section. 

Thus,  from  the  Kelvin  and  Helmholtz  theorems,  each  filament  may  be 
followed  throughout  the  time  history  of  the  flow  in  a material  reference 
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frame  with  the  circulation  T of  each  filament  remaining  constant  in 


time . 


* dt  mff 


U) 


dA 


(2.3) 


Here  A is  the  cross-section  area  of  the  filament, 
the  vorticity  field  in  each  box  is  taken  to  be 


“4  /-4 . f .*4  -4/  . -4  ,-4/ 

u)(r)  =J  G(r-r  ) u)  (r  ) 


dr 


In  particular. 


(2.4) 


where  G is  a filter  function  and  the  unfiltered  vorticity  is  generated 
by  the  space  curves  r^(5)  as  follows: 

2(?)  ri/ di  • (2.5) 

i 

Here  § is  a parameter  which  traces  each  filament  along  its  length  at 
any  instant  in  time.  The  summation  is  over  individual  vortex  filaments. 
The  evolution  of  each  space  curve  is  determined  from  the  filtered  velo- 
city by 

Sr  (§)  f 

— at — J G^i"^  ) “(^  ) (2*6) 


with  u determined  from  (2.2).  Using  the  same  filter  G in  (2.6)  as 
in  (2.4)  will  ensure  momentum  conservation. 


3.  COMPUTATIONAL  DESCRIPTION 

For  a variety  of  reasons,  the  vorticity  field  and  other  fields  are 
conveniently  expressed  in  Fourier  space.  Just  as  in  the  more  successful 
numerical  attacks  on  the  turbulence  problem  by  solution  of  the  Navier- 
Stokes  equation^^  The  main  reason  why  the  velocity  component 

is  recorded  in  spectral  form  is  that  calculus  (differentiation  and 
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Integration)  translates  to  algebra  (multiplication  and  division)  in 
the  spectral  domain.  Of  course,  this  Implies  periodicity  in  all 
three  dimensions.  However,  the  potential  of  getting  away  from  finite- 
difference  methods  by  means  of  spectra  cannot  be  fully  realized.  In 
principle,  local  evaluation  of  F(r)  by  summation  of  individually 
calculated  trigonometric  functions  sin(k*r)  , cos(k'r)  would 
eliminate  grids.  In  practice,  when  this  has  to  be  done  at  a large 
number  of  places,  as  in  vortex  tracing,  one  must  first  transform  the 
spectrum  F(k)  onto  some  grid  and  then  interpolate  for  arbitrary  r 
from  nearest  grid  data.  Likewise,  when  the  excitation  of  a spectrum 
by  local  sources  of  vorticity  is  evaluated,  a similar  act  of  inter- 
polation is  called  for. 

At  this  point,  before  going  into  details,  a schematic  description 
of  the  computation  is  given  below. 
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where  Y is  the  Fourier  transform  of  the  vector  potential.  The  differ- 
ent parts  of  the  scheme  will  now  be  explained,  namely  the  numerical 
modeling  of  the  vortex  filaments,  the  interpolations  that  take  place, 
the  shaping  of  the  vortices  and  the  time-stepping  procedure. 

3.1  Filament  Modeling:  In  our  model,  one  describes  each  vortex  fila- 
ment by  a succession  of  closely  spaced  markers.  Considering  a single 
vortex  in  (2.^),  we  have 


u!fr)  = T 6(r-r(5))  dr  <*§ 

a? 


at  an  instant  t , where  T is  given  by  (2.3).  Taking  the  Fourier 
transform,  we  get 

ttiOt)  =f  e lk  r T f 6(^-r(5))  dr  d§  dir 
If  we  now  discretize  r into  piece-wise  linear  sections, 


sot) 


(3.1) 


where  m is  the  total  number  of  markers  describing  the  filament; 

r = r . Integrating  (3.1)  and  leiting  £•  (rl'r 1-l)  = e.  . we  obtain: 
m o ~2  — J 


3(C)  - r f,  (3,-?^)  .'"“(Wn/a  ^ . 


Rather  than  interpolating  each  trigonometric  function  in  the  Fourier 
series  separately,  it  is  more  efficient  to  first  distribute  vortlcity 
onto  the  grid  according  to  an  interpolation  process  and  then  perform 


an  FFT.  To  do  this,  we  must  be  able  to  express  sin  e /e  as  com~ 
-i£-f 

bination  of  e where  f is  a function  of  r, 


■j  <md  rj-i 


r -i  o*1 

In  our  case,  evaluating  (3«l)  by  gaussian  quadrature1-  , we  obtain 


m 


(k)  = (?  -*  ) 

j=i  ? 


,-£*%  (1+3'^)^  + (l-3"^)r + e-i£,J5 


Indeed  this  is  equivalent  to  the  approximation. 


sin  £ 


-ie  3"^  ie.3"^ 

e 3 + e 3 


i 


ik*r 


What  we  want  now  is  to  replace  the  pure  harmonic  e for  arbitrary 

r by  an  approximant  to  be  evaluated  on  the  discrete  spatial  mesh. 

3.?  interpolation:  In  each  of  the  three  dimensions  we  use  quadratic 

In 

particular  e is  represented  in  the  interval  n-^^x^n+% 

as  the  superposition  of  three  parabolic  arcs  as  follows: 


ilex  ilen 

spline  interpolation  to  approximate  e in  terms  of  e 
ikx 


ikx 


, ,r.l  , , 1x2  ik  (n+1 ) J2l  ! ikn 

~ S(k)|j  (x  - n + -)  e +(£  ‘ (x  - n)  Je 


. 1 / . 3 x2  ik 

+ g(n  + 2"x)  e 


(n-3  )j 


Rather  than  choosing  the  function  S(k)  to  force  agreement  at  the 
point  x=n  , we  select  S(k)  to  minimize  the  mean  square  error  over 

T 12 1 

the  Interval  n-ifSx^n  + Vj  . This  leads  to: 


SOO-gstn  |)3/ (l  - .l.2  | + fj  .In'*  |) 

Therefore,  the  three  nearest  grid  points  in  each  dimension  (27 

in  all)  will  share  the  vorticity  distribution  according  to  the  spline 

function  weighting  on  them.  In  other  words,  this  is  a particular  way 

to  distribute  a finite-sized  vorticity  on  its  neighboring  grid  points. 

The  quadratic  spline  weighting  is  superior  to  the  zero-order  weighting 

(NGP  model)  and  first-order  weighting  (CIC  model)  in  the  sense  of 

creating  less  field-noise  and  resulting  in  smoother  simulation  func- 

[19] 

tions  . This  is  an  obvious  conclusion  since  vorticity  is  now  dis- 
tributed among  three  grid  points  instead  of  one  or  two  as  in  the  other 
models  and  the  interpolated  distribution  is  quadratic  rather  than  a 
piece-wise  step  function,  or  first-order  linear  function,  with  discon- 
tinuous derivatives.  There  is  also  a reduction  of  aliassing. 

Shaping:  It  is  intuitively  obvious  that  low-|k|  harmonics  are 

interpolated  by  a certain  tabulation  mesh  better  than  h igh  — J Vc] 
harmonics.  Aliassing  sets  a limit  at  ^max  = nM  for  each  component 
of  1c  (A  = mesh  spacing):  any  harmonic  with  a k-component  higher 
than  this  will  be  misinterpreted  by  the  interpolator  as  a corresponding 
lower  harmonic  with  all  k-components  lying  within  the  interval 
(-n/A  , tt/A)  . 

The  effects  of  a finite  cut-off  of  the  spectrum  on  the  p .ysics  to 
be  computed  is  an  important  question  separate  from  the  question  of 
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Interpolation.  Two  aspects  of  the  cut-off  problem  are  worth  emphasiz- 
ing. 

Firstly,  a sharp  cut-off  in  k-space  is  undesirable  (no  matter  how 

perfectly  each  harmonic  is  evaluated)  because  it  surrounds  the  objects 

that  interact  via  the  field  with  halos  in  r-space.  The  halos  decay 

only  weakly  with  distance,  like  1/r.  If  instead,  the  spectrum  is 

brought  to  zero  more  smoothly,  say  at  least  parabolically , such  halos 

become  attenuated  more  strongly.  A bell-shaped  cut-off  factor  G‘"(|k|) 

suggests  itself.  Consequently,  even  if  the  object  of  interpolation 

studies  is  to  push  up  the  maximum  usable  k , it  should  concentrate  on 

performance  in  the  range  of  low  and  intermediate  |k|  ;]lc| -values  near 

k become  irrelevant, 
max 

Secondly,  any  cut-off  factor  G ( | k ] ) implies  a shape  for  the 
interacting  vortices  which  is  the  Fourier  transform  of  G(]k[):  this 
is  so  because  the  shape  enters  the  field  interaction  process  twice, 
namely  both  when  the  field  harmonics  are  excited  by  the  local  sources 
(vortices)  and  when  the  field  harmonics  react  back  on  them.  So,  any 
such  factor  in  the  spectral  domain  (introduced  primarily  for  the  purpose 
of  fitting  the  field  harmonics  into  a finite  computer)  can  be  interpreted 
physically  as  vorticity-spreading. 

There  are  good  reasons  for  introducing  shapes  of  the  interacting 
elements  even  when  no  interpolation  is  used  at  all  and  spectral  data 
are  evaluated  precisely,  without  grids  and  tabulations.  Firstly,  the 
interacting  "fluid  elements"  in  the  real  world  are  usually  much  more 
numerous  than  those  that  can  be  accommodated  in  a computer  with  its 
peripheral  storage.  Each  element  in  the  computer  stands  in  for  a swarm 
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of  real  elements.  It  should  therefore  be  given  a spread.  Secondly, 
the  binary  interaction  of  "fluid  elements"  is  unrealistically  strong 
if  they  are  treated  as  delta-functions  in  space. 

We  note  that  the  shape  factor  should  be  isotropic  in  space:  it 
knows  no  coordinate  axes.  We  choose  an  approximately  gaussian  profile, 
equivalent  to  a gaussian  shape  in  (x,y , z )-space , but  brought  to  a 
strict  zero  at  some  maximum  |k|  . De-emphasizing  the  harmonics  close 
to,  and  just  below  k will  further  reduce  aliassing. 

3.4  Solving:  If  we  now  consider  equation  (2.2)  in  Fourier  space,  the 
velocity  field  at  the  location  of  a "vortex-marker"  is  obtained  by 
weighting  the  entries  in  the  table  of  spline  amplitudes  with  the  spline 
weights.  The  latter  are  deduced  from  the  relative  position  of  a vortex 
in  its  interpolation  cell.  The  spline  amplitudes  are  obtained  from  the 
velocity  harmonics  by  first  multiplying  with  a factor  S(k^)  and  two 
similar  factors  which  have  k^  and  kz  in  place  of  k^  , then  calling 
a three-dimensional  FFT  on  the  resulting  array. 

The  same  factors  appear  again  when  the  displacement  of  the  k*"*1 
vorticity  harmonic  is  calculated.  The  interpolation  can  be  done  for  the 
sum  of  all  the  harmonics,  and  the  tables  into  which  one  interpolates  are 
then  the  FFT's  of  the  harmonics  of  vorticity,  modified  by  the  factors 
insuring  best  mean  square  fit  in  each  dimension. 

In  going  from  the  table  of  spline  amplitudes  for  vorticity  to  the 
table  of  spline  amplitudes  for  velocity  harmonics,  one  therefore  has 
not  only  to  perform  a forward  and  backward  FFT,  with  equation  (2.2)  in 
k-space  ir.  between,  but  one  must  also  introduce  the  squares  of  the 
spline  fitting  factors  indicated  above. 
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Similarly,  we  mentioned  that  any  vorticlty  shaping  factor  should 
be  Introduced  both  when  the  local  velocity  field  action  on  the  distri- 
buted vorticlty  cloud  Is  evaluated  and  when  its  excitation  of  the 
vorticlty  harmonics  Is  accumulated.  In  both  cases,  one  could  perform 
a convolution  in  (x,y,z)-space,  but  is  is  much  quicker  to  replace  this 
by  a multiplication  in  k-space.  The  transform  of  the  shape  factor  is 
therefore  introduced  squared  along  with  the  above  mentioned  spline 
fitting  factors  in  the  course  of  solving  the  equations  for  the  velocity 
field  in  lt-space.  It  is  convenient  to  introduce  the  (squared)  shape 

p 

factor  along  with  the  inverse  Poisson  operator  1/k ‘ . In  our  present 

o 2 

16  code  there  are  only  64  possible  different  values  of  k in  the 
sphere  |k|  < ^max  » so  any  functlon  of  |lcj  can  readily  be  pre tabu- 
lated. 

8.5  Time-Stepping:  So  far  in  the  code,  we  have  been  using  the  well 
known  leap-frog  method  with  the  first  step  generated  by  Euler's  method. 
This  method  is  unstable^^  but  was  used  mainly  to  test  our  code.  In 
fact,  we  already  noticed  the  instability  of  this  scheme  in  our  two-ring 
experiments . 

Nevertheless,  used  with  an  occasional  forward  Euler  step,  it  seems 
to  be  possible  to  suppress  the  weak  instability  associated  with  leap- 
frog differencing^^. 

4.  RESULTS  OF  THE  COMPUTER  EXPERIMENTS 

The  computer  used  was  a CDC  76OO  located  at  NASA-Ames  Research 
Center.  The  present  mesh  size  is  l6  . 

A first  test  was  done  on  a single  vortex  ring  of  radius  R about 
the  z-axis.  Its  center  is  initially  located  at  (8,8,8)  in  our  mesh 
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and  thereafter  moves  along  the  z-axis.  The  circulation  Is  T = 2 . 

In  particular  we  investigated  the  initial  speed  of  the  vortex  ring  as 
a function  of  radius  and  position  around  the  ring. 

To  check  the  accuracy  of  our  mesh  technique  we  also  computed  the 
speed  of  the  vortex  ring  using  a continuum  or  Green's  function  approach. 
Since  the  filter  we  use  in  the  mesh  method  is  approximately  gaussian, 
we  consider  a single  vortex  ring  of  gaussian  cross-section.  Following 
the  procedure  as  defined  by  (2.4-2. 6),  we  obtain  ns  its  filtered  velo- 
city of  translation  in  free  space 


*1 

dt 


r /iiaL_ 
W M3 


(r-r  ; ) X hr  ' 

a?  5 


(4.1) 


where  a = J^-r/  | /2^cr  , o is  the  radius  of  the  cross-sect  ion  or 

. _2 

the  width  of  the  gaussian  filter  and  f(a)  = 2n  ?ae  - erf (a)  . 

P p 

(it  should  be  mentioned  that  for  o « R'  , (4.1  ) can  be  approxi- 
mated to  yield  8r  w[ P[ n( c]  ^ where  C = 1 .058  and  e?  is 

ftt  4nR  L \a  /'  J Z 

the  unit  vector  in  the  direction  of  translation  z . Note  that 
the  actual  speed  of  a thin  vortex  ring  with  a gaussian  distribution 
of  vorticity  has  been  calculated  by  SaffmanL  and  is  given  by  the 
above  formula  but  with  C = 0.558.  The  difference  is  due  to  the  fact 
that  Saffman's  result  is  based  on  the  collective  motion  of  an  infinite 
number  of  vortex  tubes  with  internal  interaction  between  the  filaments 
whereas  our  result  represents  the  speed  of  a single  computational 
ring  filament).  To  (4.1),  we  then  add  the  Riot-Savart  contributions 
of  the  periodic  images. 

In  Figures  1 and  2,  we  plot  the  total  velocity  of  translation  versus 


R for  a periodic  array  of  our  single  rings.  Our  vortex- In-cell 
results  are  compared  with  the  Green's  function  method  given  by  (4.1 ) 
plus  image  contributions.  The  gausslan  width  used  in  the  Green's 
function  calculations  was  chosen  to  give  the  best  fit  to  the  vortex- 

t 2 

ln-cell  results  and  was  found  to  be  o « 1.1  times  the  cell  area. 

2 2 

Ibis  Is  In  good  agreement  with  a theoretical  estimate  of  a = 12/n 
1.2  based  on  a gausslan  fit  to  the  low  |l?|  behavior  of  our  filter. 

Recall  that  our  filter  Is  not  strictly  gausslan  but  Is  brought  smoothly 
to  zero  at  |lc|  - n . Figure  1 shows  the  velocity  recorded  at  points 
of  the  ring  close  to  the  x-  and  y-axes,  where  the  velocity  should 
be  minimum  since  the  Images  are  closer.  Figure  2 shows  the  velocity 
recorded  at  points  45  degrees  off  the  x-  and  y-axes,  where  the 
velocity  should  be  a maximum. 

The  four  next  figures  show  pictures  of  the  initial  velocity  field 
in  the  middle  of  the  mesh  cells  for  a ring  of  radius  R = 4. 

Figures  3 and  4 show  the  field  in  the  planes  1.5  mesh  units  below 
and  above  the  plane  of  the  ring,  respectively  at  z = 6.5  and  z = 9*5. 

where  the  magnitude  of  the  field  is  the  same  but  pointing  in  opposite 

directions,  respectively  toward  and  away  from  the  center  of  the  ring. 

Figures  5 and  6 show  the  field  in  the  planes  x *=  6.5  and  x = 9*5  • 

i 

Figure  7 shows  the  lateral  vortex  profile  in  the  (x-z)-plane  at  four 
instants.  We  can  see  the  constancy  of  the  motion. 

From  these  results,  we  derived  an  estimate  of  the  CPU  time  per 

\ ] 

time  step  to  move  a vortex  made  of  m markers.  For  m 360,  it  takes 

0.4 1 CPU  seconds  per  time  step;  for  m = 7?0,  0.48  CPU  seconds  per 

time  step.  So  in  general,  it  takes  0.34  + 8econt,s  Ppr  1 

step.  All  calculations  were  done  with  leap-frog  stepping  in  time 


Figure  3 


Projection  on  (x-y)-plane  (z  « 6.5)  of  the  velocity  field 
in  the  middle  of  the  cells  for  a single  vortex  ring. 
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Figure  4.  Projection  on  (x-y)-plane  (z  * 9*5)  of  the  velocity  field 
in  the  middle  of  the  cells  for  a single  vortex  ring. 
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Figure  7.  Displacement  of  e single  vortex  ring. 


20 


that  involves  only  one  evaluation  of  the  derivative  per  step. 

A second  test  was  done  on  a set  of  two  vortex  rings  of  radius 
R = 4 about  the  z-axis.  Their  centers  are  initially  located  at 
(8,8,7)  and  (8,8,10)  in  the  mesh  and  thereafter  move  also  along  the 
z-axis.  Both  have  the  same  circulation  T = 2 . 

We  know  that  two  similar  vortex  rings  at  some  distance  apart  on 
a common  axis  of  symmetry  play  the  following  game.  The  velocity  field 
associated  with  the  rear  vortex  ring  has  a radially  outward  component 
at  the  position  of  the  front  ring  and  so  the  radius  of  the  front  ring 
gradually  increases  (with  T constant).  This  leads  to  a decrease  in 
its  velocity  of  translation,  and  there  is  a corresponding  increase  in 
the  velocity  of  translation  of  the  rear  vortex  which  ultimately  passes 
through  the  larger  vortex  and  in  turn  becomes  the  front  vortex.  The 
maneuver  is  then  repeated.  Indeed,  we  observed  that  maneuver. 

Figures  8,  9 and  10  show  the  initial  velocity  field  respectively 
in  the  planes  z = 5*5,  z = 8.5  and  z = 11.5  • As  expected,  at 
z = 5*5  and  z = 11.5  » the  magnitude  of  the  field  is  the  same  but 
pointing  in  opposite  directions,  respectively  toward  and  away  from 
the  center.  At  z = 8.5  , centrally  between  the  two  rings,  the  field 
reaches  a minimum.  Figures  11  and  12  were  taken  at  x = 6.5  and 
x = 9*5  • The  last  five  figures  (13-17)  show  the  displacement  in  the 
(x-z)-plane  and  the  (x-y)-plane  at  ten  instants.  We  see  the  rings 
going  through  each  other  repeatedly  and  the  buildup  of  distortions  due 
to  the  influence  of  images. 


21 


/ / 


s x \ \ 1 I / / 

' N \ \ \ / / / 

N \ \ \ \ / / / 

^ \ \ \ \ / / / 


///  / \ \\ 
/ / / / \ \ \ 
/ / / I \ \ \ 

///  I \ \ \ 


/ ' 


\ \ 
\ \ 


Figure  8.  Projection  on  (x-y)-plane  (z  = 5.5)  of  the  velocity  field 
in  the  middle  of  the  cells  for  a pair  of  vortex  rings. 
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Figure  31.  Projection  on  (y-z)-plane  (x  = 6.5)  of  the  velocity  field 
in  the  middle  of  the  cells  for  a pair  of  vortex  rings. 
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Figure  12.  Projection  on  (y-z)-plane  (x  = 9.5>  of  the  velocity  field 
in  the  middle  of  the  cells  for  a pair  of  vortex  rings. 


5.  FURTHER  WORK 


As  our  next  step,  we  shall  initialize  our  code  to  the  Taylor- 

[221 

Green  vortex  system  . This  system  has  the  three-dimensional 
periodic  structure  which  is  built  into  our  first  version  of  the  code 
(by  virtue  of  using  complex  FFT's).  The  Taylor-Green  system  has 

[13  <?3l 

undergone  both  analytical  study  ’ and  numerical  study  by  methods 

r 2k  25I 

other  than  oursL  ' . The  system  is  one  of  continuous  vorticity, 

and  it  is  important  to  represent  such  a continuum  by  a sufficient 

number  of  discrete  vortex  filaments.  Sensitivity  to  the  coarseness 

of  this  discretization  (which  is  distinct  from  the  discretization  of 

the  spatial  mesh,  and  the  time  stepping)  is  to  be  explored. 

One  essential  improvement  of  our  code  will  be  to  proceed  to 
o 

finer  meshes  than  16  . With  a finer  mesh  one  could  then  explore  the 
evolution  of  vortex  structures  which  have  not  been  studied  hitherto, 
and  make  reliable  predictions  from  the  computer  output. 

Several  further  code  modifications  are  called  for.  The  periodic 
boundaries  ought  to  be  changed  so  that  one  simulates  (possibly  vortex- 
shedding)  planar  walls  in  one  or  two  of  the  three  dimensions.  This 
Is  done  by  replacing  complex  FFT's  with  readily  available  sine  and 
cosine  transforms.  One  would  retain  periodicity  in  the  third  dimen- 
sion, having  in  mind  the  simulation  of  channel  flow  or  flow  through 
a re-entrant  wind-tunnel  of  rectangular  cross-section  (but  with  curva- 
ture effects  absent).  Another  variant  is  the  simulation  of  "infinite" 


boundary  conditions  using  similar  approaches  as  in 


[26] 


and  in 


[27] 
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